home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / SoundAndMusic / cmix / filters / pm.f < prev    next >
Text File  |  1988-10-03  |  21KB  |  749 lines
  1. c
  2. c-----------------------------------------------------------------------
  3. c main program: fir linear phase filter design program
  4. c
  5. c authors: james h. mcclellan
  6. c       schlumberger well services
  7. c          12125 technology blvd.
  8. c          austin, texas 78759
  9. c
  10. c          thomas w. parks
  11. c          department of electrical and computer engineering
  12. c          rice university
  13. c          houston, texas 77251-1892
  14. c
  15. c          lawrence r. rabiner
  16. c          bell laboratories
  17. c          murray hill, new jersey 07974
  18. c
  19. c modified for terminal input by t.w. parks-9/85
  20. c the filter specifications are written to the file "pmfil.lst"
  21. c the impulse response is written to the file "response.dat"
  22. c
  23. c input:
  24. c  nfilt-- filter length
  25. c  jtype-- type of filter
  26. c          1 = multiple passband/stopband filter
  27. c          2 = differentiator
  28. c          3 = hilbert transform filter
  29. c  nbands-- number of bands
  30. c  lgrid-- grid density, will be set to 16 unless
  31. c          specified otherwise by a positive constant.
  32. c
  33. c  edge(2*nbands)-- bandedge array, lower and upper edges for each band
  34. c                   with a maximum of 10 bands.
  35. c
  36. c  fx(nbands)-- desired function array (or desired slope if a
  37. c               differentiator) for each band.
  38. c
  39. c  wtx(nbands)-- weight function array in each band.  for a
  40. c                differentiator, the weight function is inversely
  41. c                proportional to f.
  42. c
  43. c  sample input data setup:
  44. c       32,1,3,0
  45. c       0.0,0.1,0.2,0.35,0.425,0.5
  46. c       0.0,1.0,0.0
  47. c       10.0,1.0,10.0
  48. c     this data specifies a length 32 bandpass filter with
  49. c     stopbands 0 to 0.1 and 0.425 to 0.5, and passband from
  50. c     0.2 to 0.35 with weighting of 10 in the stopbands and 1
  51. c     in the passband.  the grid density defaults to 16.
  52. c
  53. c     the following input data specifies a length 32 fullband
  54. c     differentiator with slope 1 and weighting of 1/f.
  55. c     the grid density will be set to 20.
  56. c       32,2,1,20
  57. c       0,0.5
  58. c       1.0
  59. c       1.0
  60. c
  61. c-----------------------------------------------------------------------
  62. c
  63.       common pi2,ad,dev,x,y,grid,des,wt,alpha,iext,nfcns,ngrid
  64.       common /oops/niter,iout
  65.       dimension iext(66),ad(66),alpha(66),x(66),y(66)
  66.       dimension h(66),hh(132)
  67.       dimension des(1045),grid(1045),wt(1045)
  68.       dimension edge(20),fx(10),wtx(10),deviat(10)
  69.       double precision pi2,pi
  70.       double precision ad,dev,x,y
  71.       double precision gee,d
  72.       integer bd1,bd2,bd3,bd4
  73.       data bd1,bd2,bd3,bd4/1hb,1ha,1hn,1hd/
  74.       input= 5
  75.       iout= 3
  76.       pi=4.0*datan(1.0d0)
  77.       pi2=2.0d00*pi
  78. c
  79. c  the program is set up for a maximum length of 128, but
  80. c  this upper limit can be changed by redimensioning the
  81. c  arrays iext, ad, alpha, x, y, h to be nfmax/2 + 2.
  82. c  the arrays des, grid, and wt must dimensioned
  83. c  16(nfmax/2 + 2).
  84. c
  85.       nfmax=128
  86.   100 continue
  87.       open(3,file='pmfil.lst')
  88.       open(4,file='response.dat')
  89.       open(18,file='cmixpm')
  90.       jtype=0
  91. c
  92. c  program input section
  93. c
  94.       write(*,104)
  95.   104 format(3x,'Enter filter length, type, no. of bands')
  96.       read(*,*) nfilt,jtype,nbands
  97.       if(nfilt.eq.0)stop
  98.       if(nfilt.le.nfmax.or.nfilt.ge.3) go to 115
  99.       call error
  100.       stop
  101.   115 if(nbands.le.0) nbands=1
  102. c
  103. c  grid density is assumed to be 16 unless specified
  104. c  otherwise
  105. c
  106.       if(lgrid.le.0) lgrid=16
  107.       jb=2*nbands
  108.       write(*,120)
  109.   120 format(3x,'Enter band edges')    
  110.       read(*,*) (edge(j),j=1,jb)
  111.       write(*,121)
  112.   121 format(3x,'Enter desired value in each band')
  113.       read (*,*) (fx(j),j=1,nbands)
  114.       write(*,122)
  115.   122 format(3x,'Enter weight in each band')
  116.       read(*,*) (wtx(j),j=1,nbands)
  117.       if(jtype.gt.0.and.jtype.le.3) go to 125
  118.       call error
  119.       stop
  120.   125 neg=1
  121.       if(jtype.eq.1) neg=0
  122.       nodd=nfilt/2
  123.       nodd=nfilt-2*nodd
  124.       nfcns=nfilt/2
  125.       if(nodd.eq.1.and.neg.eq.0) nfcns=nfcns+1
  126. c
  127. c  set up the dense grid.  the number of points in the grid
  128. c  is (filter length + 1)*grid density/2
  129. c
  130.       grid(1)=edge(1)
  131.       delf=lgrid*nfcns
  132.       delf=0.5/delf
  133.       if(neg.eq.0) go to 135
  134.       if(edge(1).lt.delf) grid(1)=delf
  135.   135 continue
  136.       j=1
  137.       l=1
  138.       lband=1
  139.   140 fup=edge(l+1)
  140.   145 temp=grid(j)
  141. c
  142. c  calculate the desired magnitude response and the weight
  143. c  function on the grid
  144. c
  145.       des(j)=eff(temp,fx,wtx,lband,jtype)
  146.       wt(j)=wate(temp,fx,wtx,lband,jtype)
  147.       j=j+1
  148.       grid(j)=temp+delf
  149.       if(grid(j).gt.fup) go to 150
  150.       go to 145
  151.   150 grid(j-1)=fup
  152.       des(j-1)=eff(fup,fx,wtx,lband,jtype)
  153.       wt(j-1)=wate(fup,fx,wtx,lband,jtype)
  154.       lband=lband+1
  155.       l=l+2
  156.       if(lband.gt.nbands) go to 160
  157.       grid(j)=edge(l)
  158.       go to 140
  159.   160 ngrid=j-1
  160.       if(neg.ne.nodd) go to 165
  161.       if(grid(ngrid).gt.(0.5-delf)) ngrid=ngrid-1
  162.   165 continue
  163. c
  164. c  set up a new approximation problem which is equivalent
  165. c  to the original problem
  166. c
  167.       if(neg) 170,170,180
  168.   170 if(nodd.eq.1) go to 200
  169.       do 175 j=1,ngrid
  170.       change=dcos(pi*grid(j))
  171.       des(j)=des(j)/change
  172.   175 wt(j)=wt(j)*change
  173.       go to 200
  174.   180 if(nodd.eq.1) go to 190
  175.       do 185 j=1,ngrid
  176.       change=dsin(pi*grid(j))
  177.       des(j)=des(j)/change
  178.   185 wt(j)=wt(j)*change
  179.       go to 200
  180.   190 do 195 j=1,ngrid
  181.       change=dsin(pi2*grid(j))
  182.       des(j)=des(j)/change
  183.   195 wt(j)=wt(j)*change
  184. c
  185. c  initial guess for the extremal frequencies--equally
  186. c  spaced along the grid
  187. c
  188.   200 temp=float(ngrid-1)/float(nfcns)
  189.       do 210 j=1,nfcns
  190.       xt=j-1
  191.   210 iext(j)=xt*temp+1.0
  192.       iext(nfcns+1)=ngrid
  193.       nm1=nfcns-1
  194.       nz=nfcns+1
  195. c
  196. c  call the remes exchange algorithm to do the approximation
  197. c  problem
  198. c
  199.       call remes
  200. c
  201. c  calculate the impulse response.
  202. c
  203.       if(neg) 300,300,320
  204.   300 if(nodd.eq.0) go to 310
  205.       do 305 j=1,nm1
  206.       nzmj=nz-j
  207.   305 h(j)=0.5*alpha(nzmj)
  208.       h(nfcns)=alpha(1)
  209.       go to 350
  210.   310 h(1)=0.25*alpha(nfcns)
  211.       do 315 j=2,nm1
  212.       nzmj=nz-j
  213.       nf2j=nfcns+2-j
  214.   315 h(j)=0.25*(alpha(nzmj)+alpha(nf2j))
  215.       h(nfcns)=0.5*alpha(1)+0.25*alpha(2)
  216.       go to 350
  217.   320 if(nodd.eq.0) go to 330
  218.       h(1)=0.25*alpha(nfcns)
  219.       h(2)=0.25*alpha(nm1)
  220.       do 325 j=3,nm1
  221.       nzmj=nz-j
  222.       nf3j=nfcns+3-j
  223.   325 h(j)=0.25*(alpha(nzmj)-alpha(nf3j))
  224.       h(nfcns)=0.5*alpha(1)-0.25*alpha(3)
  225.       h(nz)=0.0
  226.       go to 350
  227.   330 h(1)=0.25*alpha(nfcns)
  228.       do 335 j=2,nm1
  229.       nzmj=nz-j
  230.       nf2j=nfcns+2-j
  231.   335 h(j)=0.25*(alpha(nzmj)-alpha(nf2j))
  232.       h(nfcns)=0.5*alpha(1)-0.25*alpha(2)
  233. c
  234. c  program output section.
  235. c
  236. c since iout=3, the output is written to file 'pmfil.lst'
  237.   350 write(iout,360)
  238.   360 format(1h1, 70(1h*)/15x,29hfinite impulse response (fir)/
  239.      113x,34hlinear phase digital filter design/
  240.      217x,24hremes exchange algorithm)
  241.       if(jtype.eq.1) write(iout,365)
  242.   365 format(22x,15hbandpass filter/)
  243.       if(jtype.eq.2) write(iout,370)
  244.   370 format(22x,14hdifferentiator/)
  245.       if(jtype.eq.3) write(iout,375)
  246.   375 format(20x,19hhilbert transformer/)
  247.       write(iout,378) nfilt
  248.   378 format(20x,16hfilter length = ,i3/)
  249. c for screen output, the impulse response is not written
  250.       if (iout.eq.6) go to 457
  251.       write(iout,380)
  252.   380 format(15x,28h***** impulse response *****)
  253.       do 381 j=1,nfcns
  254.       k=nfilt+1-j
  255.       if(neg.eq.0) write(iout,382) j,h(j),k
  256.       if(neg.eq.1) write(iout,383) j,h(j),k
  257.   381 continue
  258. c now to write impulse response to file 'response.dat'
  259. c write the first half of the response
  260. c********************* write no. of coeffs. for convenience *******
  261.       write(4,7865)nfilt
  262.  7865 format(i5)
  263.       do 785 i=1,nfcns
  264.       hh(i)=h(i)
  265.   785 continue
  266.       if(nodd.lt.1) then
  267.       do 786 n=1,nfcns
  268.       hh(nfcns+n)=h(nfcns-n+1)
  269.   786 continue
  270.       do 787 i=1,nfcns*2
  271.       write(4,783) hh(i)
  272.   787 continue
  273.       call cmix(hh,nfcns*2)
  274.   783 format(g20.5)
  275.       else
  276.       do 788 n=1,nfcns-1 
  277.       hh(nfcns+n)=h(nfcns-n)
  278.   788 continue
  279.       do 789 i=1,nfcns*2-1
  280.       write(4,783) hh(i)
  281.   789 continue
  282.       call cmix(hh,nfcns*2-1)
  283.       end if
  284.   382 format(13x,2hh(,i2,4h) = ,e15.8,5h = h(,i3,1h))
  285.   383 format(13x,2hh(,i2,4h) = ,e15.8,6h = -h(,i3,1h))
  286.       if(neg.eq.1.and.nodd.eq.1) write(iout,384) nz
  287.   384 format(13x,2hh(,i2,8h) =  0.0)
  288.   457 do 450 k=1,nbands,4
  289.       kup=k+3
  290.       if(kup.gt.nbands) kup=nbands
  291.       write(iout,385) (bd1,bd2,bd3,bd4,j,j=k,kup)
  292.   385 format(24x,4(4a1,i3,7x))
  293.       write(iout,390) (edge(2*j-1),j=k,kup)
  294.   390 format(2x,15hlower band edge,5f14.7)
  295.       write(iout,395) (edge(2*j),j=k,kup)
  296.   395 format(2x,15hupper band edge,5f14.7)
  297.       if(jtype.ne.2) write(iout,400) (fx(j),j=k,kup)
  298.   400 format(2x,13hdesired value,2x,5f14.7)
  299.       if(jtype.eq.2) write(iout,405) (fx(j),j=k,kup)
  300.   405 format(2x,13hdesired slope,2x,5f14.7)
  301.       write(iout,410) (wtx(j),j=k,kup)
  302.   410 format(2x,9hweighting,6x,5f14.7)
  303.       do 420 j=k,kup
  304.   420 deviat(j)=dev/wtx(j)
  305.       write(iout,425) (deviat(j),j=k,kup)
  306.   425 format(2x,9hdeviation,6x,5f14.7)
  307.       if(jtype.ne.1) go to 450
  308.       do 430 j=k,kup
  309.   430 deviat(j)=20.0*alog10(abs(deviat(j)+fx(j)))
  310.       write(iout,435) (deviat(j),j=k,kup)
  311.   435 format(2x,15hdeviation in db,5f14.7)
  312.   450 continue
  313.       do 452 j=1,nz
  314.       ix=iext(j)
  315.   452 grid(j)=grid(ix)
  316. c extremal frequencies not written out in this version
  317. c      write(iout,455) (grid(j),j=1,nz)
  318.   455 format(/2x,47hextremal frequencies--maxima of the error curve/
  319.      1 (2x,5f12.7))
  320. c      write(iout,460)
  321.   460 format(1x,70(1h*)/)
  322.       if (iout.eq.6) go to 461
  323.       iout=6
  324.       go to 350
  325.   461 write(*,460)
  326.       write(*,462)
  327.   462 format(10x,48h*****filter specs are in the file pmfil.lst*****/)
  328.   458 write(*,459)
  329.   459 format(10x,48h****impulse response is in file response.dat****/)
  330.       write(*,463)
  331.   463 format(10x,'***cmix data is in file cmixpm****')
  332.       stop
  333.       end
  334. c
  335. c-----------------------------------------------------------------------
  336. c function: eff
  337. c   function to calculate the desired magnitude response
  338. c   as a function of frequency.
  339. c   an arbitrary function of frequency can be
  340. c   approximated if the user replaces this function
  341. c   with the appropriate code to evaluate the ideal
  342. c   magnitude.  note that the parameter freq is the
  343. c   value of normalized frequency needed for evaluation.
  344. c-----------------------------------------------------------------------
  345. c
  346.       function eff(freq,fx,wtx,lband,jtype)
  347.       dimension fx(5),wtx(5)
  348.       if(jtype.eq.2) go to 1
  349.       eff=fx(lband)
  350.       return
  351.     1 eff=fx(lband)*freq
  352.       return
  353.       end
  354. c
  355. c-----------------------------------------------------------------------
  356. c function: wate
  357. c   function to calculate the weight function as a function
  358. c   of frequency.  similar to the function eff, this function can
  359. c   be replaced by a user-written routine to calculate any
  360. c   desired weighting function.
  361. c-----------------------------------------------------------------------
  362. c
  363.       function wate(freq,fx,wtx,lband,jtype)
  364.       dimension fx(5),wtx(5)
  365.       if(jtype.eq.2) go to 1
  366.       wate=wtx(lband)
  367.       return
  368.     1 if(fx(lband).lt.0.0001) go to 2
  369.       wate=wtx(lband)/freq
  370.       return
  371.     2 wate=wtx(lband)
  372.       return
  373.       end
  374. c
  375. c-----------------------------------------------------------------------
  376. c subroutine: error
  377. c   this routine writes an error message if an
  378. c   error has been detected in the input data.
  379. c-----------------------------------------------------------------------
  380. c
  381.       subroutine error
  382.       common /oops/niter,iout
  383.       write(iout,1)
  384.     1 format(44h ************ error in input data **********)
  385.       return
  386.       end
  387. c
  388. c-----------------------------------------------------------------------
  389. c subroutine: remes
  390. c   this subroutine implements the remes exchange algorithm
  391. c   for the weighted chebyshev approximation of a continuous
  392. c   function with a sum of cosines.  inputs to the subroutine
  393. c   are a dense grid which replaces the frequency axis, the
  394. c   desired function on this grid, the weight function on the
  395. c   grid, the number of cosines, and an initial guess of the
  396. c   extremal frequencies.  the program minimizes the chebyshev
  397. c   error by determining the best location of the extremal
  398. c   frequencies (points of maximum error) and then calculates
  399. c   the coefficients of the best approximation.
  400. c-----------------------------------------------------------------------
  401. c
  402.       subroutine remes
  403.       common pi2,ad,dev,x,y,grid,des,wt,alpha,iext,nfcns,ngrid
  404.       common /oops/niter,iout
  405.       dimension iext(66),ad(66),alpha(66),x(66),y(66)
  406.       dimension des(1045),grid(1045),wt(1045)
  407.       dimension a(66),p(65),q(65)
  408.       double precision pi2,dnum,dden,dtemp,a,p,q
  409.       double precision dk,dak
  410.       double precision ad,dev,x,y
  411.       double precision gee,d
  412. c
  413. c  the program allows a maximum number of iterations of 25
  414. c
  415.       itrmax=25
  416.       devl=-1.0
  417.       nz=nfcns+1
  418.       nzz=nfcns+2
  419.       niter=0
  420.   100 continue
  421.       iext(nzz)=ngrid+1
  422.       niter=niter+1
  423.       if(niter.gt.itrmax) go to 400
  424.       do 110 j=1,nz
  425.       jxt=iext(j)
  426.       dtemp=grid(jxt)
  427.       dtemp=dcos(dtemp*pi2)
  428.   110 x(j)=dtemp
  429.       jet=(nfcns-1)/15+1
  430.       do 120 j=1,nz
  431.   120 ad(j)=d(j,nz,jet)
  432.       dnum=0.0
  433.       dden=0.0
  434.       k=1
  435.       do 130 j=1,nz
  436.       l=iext(j)
  437.       dtemp=ad(j)*des(l)
  438.       dnum=dnum+dtemp
  439.       dtemp=float(k)*ad(j)/wt(l)
  440.       dden=dden+dtemp
  441.   130 k=-k
  442.       dev=dnum/dden
  443. c      write(iout,131) dev
  444. c intermeditate deviations not written in this version
  445.   131 format(1x,12hdeviation = ,f12.9)
  446.       nu=1
  447.       if(dev.gt.0.0) nu=-1
  448.       dev=-float(nu)*dev
  449.       k=nu
  450.       do 140 j=1,nz
  451.       l=iext(j)
  452.       dtemp=float(k)*dev/wt(l)
  453.       y(j)=des(l)+dtemp
  454.   140 k=-k
  455.       if(dev.gt.devl) go to 150
  456.       call ouch
  457.       go to 400
  458.   150 devl=dev
  459.       jchnge=0
  460.       k1=iext(1)
  461.       knz=iext(nz)
  462.       klow=0
  463.       nut=-nu
  464.       j=1
  465. c
  466. c  search for the extremal frequencies of the best
  467. c  approximation
  468. c
  469.   200 if(j.eq.nzz) ynz=comp
  470.       if(j.ge.nzz) go to 300
  471.       kup=iext(j+1)
  472.       l=iext(j)+1
  473.       nut=-nut
  474.       if(j.eq.2) y1=comp
  475.       comp=dev
  476.       if(l.ge.kup) go to 220
  477.       err=gee(l,nz)
  478.       err=(err-des(l))*wt(l)
  479.       dtemp=float(nut)*err-comp
  480.       if(dtemp.le.0.0) go to 220
  481.       comp=float(nut)*err
  482.   210 l=l+1
  483.       if(l.ge.kup) go to 215
  484.       err=gee(l,nz)
  485.       err=(err-des(l))*wt(l)
  486.       dtemp=float(nut)*err-comp
  487.       if(dtemp.le.0.0) go to 215
  488.       comp=float(nut)*err
  489.       go to 210
  490.   215 iext(j)=l-1
  491.       j=j+1
  492.       klow=l-1
  493.       jchnge=jchnge+1
  494.       go to 200
  495.   220 l=l-1
  496.   225 l=l-1
  497.       if(l.le.klow) go to 250
  498.       err=gee(l,nz)
  499.       err=(err-des(l))*wt(l)
  500.       dtemp=float(nut)*err-comp
  501.       if(dtemp.gt.0.0) go to 230
  502.       if(jchnge.le.0) go to 225
  503.       go to 260
  504.   230 comp=float(nut)*err
  505.   235 l=l-1
  506.       if(l.le.klow) go to 240
  507.       err=gee(l,nz)
  508.       err=(err-des(l))*wt(l)
  509.       dtemp=float(nut)*err-comp
  510.       if(dtemp.le.0.0) go to 240
  511.       comp=float(nut)*err
  512.       go to 235
  513.   240 klow=iext(j)
  514.       iext(j)=l+1
  515.       j=j+1
  516.       jchnge=jchnge+1
  517.       go to 200
  518.   250 l=iext(j)+1
  519.       if(jchnge.gt.0) go to 215
  520.   255 l=l+1
  521.       if(l.ge.kup) go to 260
  522.       err=gee(l,nz)
  523.       err=(err-des(l))*wt(l)
  524.       dtemp=float(nut)*err-comp
  525.       if(dtemp.le.0.0) go to 255
  526.       comp=float(nut)*err
  527.       go to 210
  528.   260 klow=iext(j)
  529.       j=j+1
  530.       go to 200
  531.   300 if(j.gt.nzz) go to 320
  532.       if(k1.gt.iext(1)) k1=iext(1)
  533.       if(knz.lt.iext(nz)) knz=iext(nz)
  534.       nut1=nut
  535.       nut=-nu
  536.       l=0
  537.       kup=k1
  538.       comp=ynz*(1.00001)
  539.       luck=1
  540.   310 l=l+1
  541.       if(l.ge.kup) go to 315
  542.       err=gee(l,nz)
  543.       err=(err-des(l))*wt(l)
  544.       dtemp=float(nut)*err-comp
  545.       if(dtemp.le.0.0) go to 310
  546.       comp=float(nut)*err
  547.       j=nzz
  548.       go to 210
  549.   315 luck=6
  550.       go to 325
  551.   320 if(luck.gt.9) go to 350
  552.       if(comp.gt.y1) y1=comp
  553.       k1=iext(nzz)
  554.   325 l=ngrid+1
  555.       klow=knz
  556.       nut=-nut1
  557.       comp=y1*(1.00001)
  558.   330 l=l-1
  559.       if(l.le.klow) go to 340
  560.       err=gee(l,nz)
  561.       err=(err-des(l))*wt(l)
  562.       dtemp=float(nut)*err-comp
  563.       if(dtemp.le.0.0) go to 330
  564.       j=nzz
  565.       comp=float(nut)*err
  566.       luck=luck+10
  567.       go to 235
  568.   340 if(luck.eq.6) go to 370
  569.       do 345 j=1,nfcns
  570.       nzzmj=nzz-j
  571.       nzmj=nz-j
  572.   345 iext(nzzmj)=iext(nzmj)
  573.       iext(1)=k1
  574.       go to 100
  575.   350 kn=iext(nzz)
  576.       do 360 j=1,nfcns
  577.   360 iext(j)=iext(j+1)
  578.       iext(nz)=kn
  579.       go to 100
  580.   370 if(jchnge.gt.0) go to 100
  581. c
  582. c  calculation of the coefficients of the best approximation
  583. c  using the inverse discrete fourier transform
  584. c
  585.   400 continue
  586.       nm1=nfcns-1
  587.       fsh=1.0e-06
  588.       gtemp=grid(1)
  589.       x(nzz)=-2.0
  590.       cn=2*nfcns-1
  591.       delf=1.0/cn
  592.       l=1
  593.       kkk=0
  594.       if(grid(1).lt.0.01.and.grid(ngrid).gt.0.49) kkk=1
  595.       if(nfcns.le.3) kkk=1
  596.       if(kkk.eq.1) go to 405
  597.       dtemp=dcos(pi2*grid(1))
  598.       dnum=dcos(pi2*grid(ngrid))
  599.       aa=2.0/(dtemp-dnum)
  600.       bb=-(dtemp+dnum)/(dtemp-dnum)
  601.   405 continue
  602.       do 430 j=1,nfcns
  603.       ft=j-1
  604.       ft=ft*delf
  605.       xt=dcos(pi2*ft)
  606.       if(kkk.eq.1) go to 410
  607.       xt=(xt-bb)/aa
  608.       xt1=sqrt(1.0-xt*xt)
  609.       ft=atan2(xt1,xt)/pi2
  610.   410 xe=x(l)
  611.       if(xt.gt.xe) go to 420
  612.       if((xe-xt).lt.fsh) go to 415
  613.       l=l+1
  614.       go to 410
  615.   415 a(j)=y(l)
  616.       go to 425
  617.   420 if((xt-xe).lt.fsh) go to 415
  618.       grid(1)=ft
  619.       a(j)=gee(1,nz)
  620.   425 continue
  621.       if(l.gt.1) l=l-1
  622.   430 continue
  623.       grid(1)=gtemp
  624.       dden=pi2/cn
  625.       do 510 j=1,nfcns
  626.       dtemp=0.0
  627.       dnum=j-1
  628.       dnum=dnum*dden
  629.       if(nm1.lt.1) go to 505
  630.       do 500 k=1,nm1
  631.       dak=a(k+1)
  632.       dk=k
  633.   500 dtemp=dtemp+dak*dcos(dnum*dk)
  634.   505 dtemp=2.0*dtemp+a(1)
  635.   510 alpha(j)=dtemp
  636.       do 550 j=2,nfcns
  637.   550 alpha(j)=2.0*alpha(j)/cn
  638.       alpha(1)=alpha(1)/cn
  639.       if(kkk.eq.1) go to 545
  640.       p(1)=2.0*alpha(nfcns)*bb+alpha(nm1)
  641.       p(2)=2.0*aa*alpha(nfcns)
  642.       q(1)=alpha(nfcns-2)-alpha(nfcns)
  643.       do 540 j=2,nm1
  644.       if(j.lt.nm1) go to 515
  645.       aa=0.5*aa
  646.       bb=0.5*bb
  647.   515 continue
  648.       p(j+1)=0.0
  649.       do 520 k=1,j
  650.       a(k)=p(k)
  651.   520 p(k)=2.0*bb*a(k)
  652.       p(2)=p(2)+a(1)*2.0*aa
  653.       jm1=j-1
  654.       do 525 k=1,jm1
  655.   525 p(k)=p(k)+q(k)+aa*a(k+1)
  656.       jp1=j+1
  657.       do 530 k=3,jp1
  658.   530 p(k)=p(k)+aa*a(k-1)
  659.       if(j.eq.nm1) go to 540
  660.       do 535 k=1,j
  661.   535 q(k)=-a(k)
  662.       nf1j=nfcns-1-j
  663.       q(1)=q(1)+alpha(nf1j)
  664.   540 continue
  665.       do 543 j=1,nfcns
  666.   543 alpha(j)=p(j)
  667.   545 continue
  668.       if(nfcns.gt.3) return
  669.       alpha(nfcns+1)=0.0
  670.       alpha(nfcns+2)=0.0
  671.       return
  672.       end
  673. c
  674. c-----------------------------------------------------------------------
  675. c function: d
  676. c   function to calculate the lagrange interpolation
  677. c   coefficients for use in the function gee.
  678. c-----------------------------------------------------------------------
  679. c
  680.       double precision function d(k,n,m)
  681.       common pi2,ad,dev,x,y,grid,des,wt,alpha,iext,nfcns,ngrid
  682.       dimension iext(66),ad(66),alpha(66),x(66),y(66)
  683.       dimension des(1045),grid(1045),wt(1045)
  684.       double precision ad,dev,x,y
  685.       double precision q
  686.       double precision pi2
  687.       d=1.0
  688.       q=x(k)
  689.       do 3 l=1,m
  690.       do 2 j=l,n,m
  691.       if(j-k)1,2,1
  692.     1 d=2.0*d*(q-x(j))
  693.     2 continue
  694.     3 continue
  695.       d=1.0/d
  696.       return
  697.       end
  698. c
  699. c-----------------------------------------------------------------------
  700. c function: gee
  701. c   function to evaluate the frequency response using the
  702. c   lagrange interpolation formula in the barycentric form
  703. c-----------------------------------------------------------------------
  704. c
  705.       double precision function gee(k,n)
  706.       common pi2,ad,dev,x,y,grid,des,wt,alpha,iext,nfcns,ngrid
  707.       dimension iext(66),ad(66),alpha(66),x(66),y(66)
  708.       dimension des(1045),grid(1045),wt(1045)
  709.       double precision p,c,d,xf
  710.       double precision pi2
  711.       double precision ad,dev,x,y
  712.       p=0.0
  713.       xf=grid(k)
  714.       xf=dcos(pi2*xf)
  715.       d=0.0
  716.       do 1 j=1,n
  717.       c=xf-x(j)
  718.       c=ad(j)/c
  719.       d=d+c
  720.     1 p=p+c*y(j)
  721.       gee=p/d
  722.       return
  723.       end
  724. c
  725. c-----------------------------------------------------------------------
  726. c subroutine: ouch
  727. c   writes an error message when the algorithm fails to
  728. c   converge.  there seem to be two conditions under which
  729. c   the algorithm fails to converge: (1) the initial
  730. c   guess for the extremal frequencies is so poor that
  731. c   the exchange iteration cannot get started, or
  732. c   (2) near the termination of a correct design,
  733. c   the deviation decreases due to rounding errors
  734. c   and the program stops.  in this latter case the
  735. c   filter design is probably acceptable, but should
  736. c   be checked by computing a frequency response.
  737. c-----------------------------------------------------------------------
  738. c
  739.       subroutine ouch
  740.       common /oops/niter,iout
  741.       write(iout,1)niter
  742.     1 format(44h ************ failure to converge **********/
  743.      141h0probable cause is machine rounding error/
  744.      223h0number of iterations =,i4/
  745.      339h0if the number of iterations exceeds 3,/
  746.      462h0the design may be correct, but should be verified with an fft)
  747.       return
  748.       end
  749.